library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tibble)
library(stringr)
library(magrittr)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
##
## extract
library(ggplot2)
library(ggpubr)
library(patchwork)
Load the motor neuron counting data.
files <- list.files("~/spinal_cord_paper/data/MN_counting/", pattern = ".txt")
tables <- list()
for (i in seq(files)) {
tables[[i]] <- read.delim(paste0("~/spinal_cord_paper/data/MN_counting/", files[i]), sep = ",")
}
names(tables) <- files
Combine the tables into a single data frame and bring it into long format. Also define a split to separate sections into rostral and caudal. The value used is the median of the index range by day.
df <- do.call(rbind,tables) %>%
rownames_to_column("tmp") %>%
mutate(index = str_extract(tmp, "\\d{1,2}$")) %>%
mutate(day = str_extract(tmp, "^Day\\d{1,2}")) %>%
mutate(day = factor(day, levels = c("Day6", "Day8", "Day10"))) %>%
mutate(embryo = str_extract(tmp, "emb_[abcdef]")) %>%
mutate(index = as.integer(index)) %>%
select(ctrl, poly, index, day, embryo) %>%
mutate(embryo_id = paste(day, embryo, sep = "_"))
# factor levels to order the embryo ids
id_levels <- c("Day6_emb_a","Day6_emb_b","Day6_emb_c","Day6_emb_d","Day6_emb_e",
"Day8_emb_a","Day8_emb_b","Day8_emb_c","Day8_emb_d","Day8_emb_e",
"Day10_emb_a","Day10_emb_b","Day10_emb_c","Day10_emb_d", "Day8_emb_f")
long_df <- gather(df, "condition", "count", -index, -day, -embryo, -embryo_id) %>%
mutate(condition = factor(condition, levels = c("ctrl", "poly"))) %>%
mutate(embryo_id = factor(embryo_id, levels = id_levels))
# use median to split indeces into rostral and caudal
axis_split <- long_df %>%
group_by(day) %>%
summarise(split = median(range(index)))
splits <- axis_split$split
names(splits) <- axis_split$day
long_df <- long_df %>%
mutate(axis = case_when(
day == names(splits)[1] & index < splits[[1]] ~ "caudal",
day == names(splits)[2] & index < splits[[2]] ~ "caudal",
day == names(splits)[3] & index < splits[[3]] ~ "caudal",
TRUE ~ "rostral"
)) %>%
mutate(split = case_when(
day == names(splits)[1] ~ splits[[1]],
day == names(splits)[2] ~ splits[[2]],
day == names(splits)[3] ~ splits[[3]]
))
Plot the raw data split by embryonic days.
total counts by embryo
counts split by condition by embryo
counts split by condition
ggplot(data = long_df, aes(x = index, y = count, color = embryo)) +
geom_point() +
geom_smooth() +
geom_vline(aes(xintercept = split), lty = "dashed") +
facet_wrap("day", nrow = 3)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 52 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(data = long_df, aes(x = index, y = count, color = condition)) +
geom_point() +
geom_smooth() +
geom_vline(aes(xintercept = split), lty = "dashed") +
facet_wrap("embryo_id", nrow = 3)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 52 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(data = long_df, aes(x = index, y = count, color = condition)) +
geom_point() +
geom_smooth() +
geom_vline(aes(xintercept = split), lty = "dashed")+
facet_wrap("day", nrow = 3)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 52 rows containing missing values or values outside the scale range
## (`geom_point()`).
There is an outlier in Day10, Emb_d, slide 54 with 550 nuclei. Most likely an added 0 (typo). Since it is almost 7X the mean of values from slides 49 to 59 of all embryos and conditions (79.29487) it is removed.
Also embryo f from day 8 is removed since it was only counted partially.
mean(long_df$count[long_df$day == "Day10" &
long_df$index %in% c(49:59)],
na.rm = TRUE)
## [1] 79.29487
# set the outlier to NA
long_df$count[long_df$count == 550] <- NA
# remove Day 8 embryo f
long_df <- long_df[!(long_df$embryo_id == "Day8_emb_f"),]
ggplot(data = long_df, aes(x = index, y = count, color = embryo)) +
geom_point() +
geom_smooth() +
geom_vline(aes(xintercept = split), lty = "dashed")+
facet_wrap("day", nrow = 1)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 53 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 53 rows containing missing values or values outside the scale range
## (`geom_point()`).
sup <- ggplot(data = long_df, aes(x = index, y = count, color = condition)) +
geom_point(size = 1, alpha = 0.5) +
geom_smooth(linewidth = 0.5) +
geom_vline(aes(xintercept = split), lty = "dashed")+
scale_x_reverse() +
scale_color_manual(values = c("black", "darkgoldenrod3")) +
facet_wrap("embryo_id", nrow = 3)
main <- ggplot(data = long_df, aes(x = index, y = count, color = condition, shape = condition)) +
geom_point(size = 1, alpha = 0.2) +
geom_vline(aes(xintercept = split), lty = "dashed")+
scale_shape_manual(values = c(1, 4)) +
scale_color_manual(values = c("black", "darkgoldenrod3")) +
scale_x_reverse() +
geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs")) +
facet_wrap("day", nrow = 3) + theme_bw()
sup
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 53 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_point()`).
main
## Warning: Removed 53 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_point()`).
Next we plot barplots of total motor neuron count to see overal changes.
nMN_df <- long_df %>%
unite(embryo_id, condition, col = "id", sep = "-") %>%
group_by(id) %>%
summarize(nMN = sum(count, na.rm = TRUE)) %>%
separate(id, sep = "-", into = c("embryo_id", "cond")) %>%
mutate(day = str_extract(embryo_id, "^Day\\d{1,2}")) %>%
mutate(grouping = paste(day, cond, sep = "_")) %>%
group_by(grouping) %>%
mutate(day = factor(day, levels = c("Day6", "Day8", "Day10")))
std.error <- function(x) sd(x)/sqrt(length(x))
se_df <- nMN_df %>% summarise(se = std.error(nMN))
avg_df <- nMN_df %>% summarise(mean = mean(nMN)) %>%
left_join(se_df, by = "grouping") %>%
separate(grouping, sep = "_", into = c("day", "cond")) %>%
mutate(day = factor(day, levels = c("Day6", "Day8", "Day10")))
# Define the top and bottom of the errorbars
limits <- aes(ymax = mean + se, ymin = mean - se)
dodge <- position_dodge2(width=0.9)
# total counts ctrl vs poly by day
bar_cond <- ggplot(data = avg_df,
aes(x = day,
y = mean,
color = cond)) +
geom_bar(position = dodge, stat="identity", fill = NA) +
geom_linerange(limits, position=dodge) +
scale_color_manual(values = c("black", "darkgoldenrod3")) +
facet_wrap("day", nrow = 3, scales = "free_x") +
theme_bw()
bar_cond
compare_means(nMN ~ cond, data = nMN_df %>% ungroup() %>% filter(day == "Day6"))
compare_means(nMN ~ cond, data = nMN_df %>% ungroup() %>% filter(day == "Day8"))
compare_means(nMN ~ cond, data = nMN_df %>% ungroup() %>% filter(day == "Day10"))
compare_means(nMN ~ day,
data = nMN_df %>%
ungroup() %>%
filter(cond == "ctrl"))
my_comparisons <- list(c("Day6", "Day8"), c("Day8", "Day10"), c("Day6", "Day10"))
box_day_ctrl <- ggplot(data = nMN_df %>%
ungroup() %>%
filter(cond == "ctrl"),
aes(x = day, y = nMN, color = cond, shape = cond)) +
geom_boxplot() +
scale_color_manual(values = c("black")) +
stat_compare_means(comparisons = my_comparisons) +
geom_point() +
scale_shape_manual(values = c(1)) +
ylim(4750,9050) +
theme_bw()
box_day_poly <- ggplot(data = nMN_df %>%
ungroup() %>%
filter(cond == "poly"),
aes(x = day, y = nMN, color = cond, shape = cond)) +
geom_boxplot() +
scale_color_manual(values = c("darkgoldenrod3")) +
stat_compare_means(comparisons = my_comparisons) +
geom_point() +
scale_shape_manual(values = c(4)) +
ylim(4750,9050) +
theme_bw()
box_day <- box_day_ctrl +
box_day_poly +
plot_layout(guides = "collect") +
plot_annotation( title = "MN counts across devel")
box_day
axis_nMN <- long_df %>%
mutate(axis = case_when(
day == names(splits)[1] & index < splits[[1]] ~ "caudal",
day == names(splits)[2] & index < splits[[2]] ~ "caudal",
day == names(splits)[3] & index < splits[[3]] ~ "caudal",
TRUE ~ "rostral"
)) %>%
unite(embryo_id, condition, axis, col = "id", sep = "-") %>%
group_by(id) %>%
summarize(nMN = sum(count, na.rm = TRUE)) %>%
separate(id, sep = "-", into = c("embryo_id", "cond", "axis")) %>%
mutate(day = str_extract(embryo_id, "^Day\\d{1,2}")) %>%
mutate(grouping = paste(day, cond, axis, sep = "_")) %>%
group_by(grouping) %>%
mutate(day = factor(day, levels = c("Day6", "Day8", "Day10")))
# grouping for the barplot
axis_levels <- c("Day6_ctrl_rostral",
"Day6_poly_rostral",
"Day6_ctrl_caudal",
"Day6_poly_caudal",
"Day8_ctrl_rostral",
"Day8_poly_rostral",
"Day8_ctrl_caudal",
"Day8_poly_caudal",
"Day10_ctrl_rostral",
"Day10_poly_rostral",
"Day10_ctrl_caudal",
"Day10_poly_caudal")
se_axis_df <- axis_nMN %>% summarise(se = std.error(nMN))
avg_axis_df <- axis_nMN %>% summarise(mean = mean(nMN)) %>%
left_join(se_axis_df, by = "grouping") %>%
separate(grouping, sep = "_", into = c("day", "cond", "axis"), remove = FALSE) %>%
mutate(day = factor(day, levels = c("Day6", "Day8", "Day10"))) %>%
mutate(grouping = factor(grouping, levels = axis_levels))
# Define the top and bottom of the errorbars
limits <- aes(ymax = mean + se, ymin = mean - se)
bar_axis <- ggplot(data = avg_axis_df,
aes(x = grouping,
y = mean,
color = cond)) +
geom_bar(position = dodge, stat="identity", fill = NA) +
geom_point(data = axis_nMN,
aes(x = grouping,
y = nMN,
shape = cond)) +
geom_linerange(limits, position=dodge) +
scale_color_manual(values = c("black", "darkgoldenrod3")) +
scale_shape_manual(values = c(1, 4)) +
facet_wrap("day", nrow = 3, scales = "free_x") +
theme_bw()
bar_axis
boxpl_axis <- axis_nMN %>%
mutate(grouping = factor(grouping, levels = axis_levels)) %>%
ggplot(aes(x = grouping,
y = nMN,
color = cond)) +
geom_boxplot() +
scale_color_manual(values = c("black", "darkgoldenrod3")) +
facet_wrap("day", nrow = 3, scales = "free_x") +
theme_bw()
boxpl_axis
boxpl_0_axis <- axis_nMN %>%
mutate(grouping = factor(grouping, levels = axis_levels)) %>%
ggplot(aes(x = grouping,
y = nMN,
color = cond)) +
geom_boxplot() +
ylim(c(0,4700)) +
scale_color_manual(values = c("black", "darkgoldenrod3")) +
facet_wrap("day", nrow = 3, scales = "free_x") +
theme_bw()
boxpl_0_axis
my_comparisons <- list(c("Day6_ctrl_rostral", "Day6_poly_rostral"),
c("Day6_ctrl_caudal", "Day6_poly_caudal"),
c("Day8_ctrl_rostral", "Day8_poly_rostral"),
c("Day8_ctrl_caudal", "Day8_poly_caudal"),
c("Day10_ctrl_rostral","Day10_poly_rostral"),
c("Day10_ctrl_caudal", "Day10_poly_caudal")
)
# total counts per day by condition
box_axis <- ggplot(data = axis_nMN %>%
ungroup() %>%
mutate(grouping = factor(grouping, levels = axis_levels)),
aes(x = grouping, y = nMN, color = cond)) +
geom_boxplot() +
scale_color_manual(values = c("black", "darkgoldenrod3")) +
stat_compare_means(comparisons = my_comparisons) +
ggtitle("rostral vs caudal MN counts across devel") +
facet_wrap("day", nrow = 3) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
box_axis
ggsave(
filename = "~/spinal_cord_paper/figures/line_plot_MN_count_individual.pdf",
width = 10, height = 6,
plot = sup
)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 53 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 53 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(
filename = "~/spinal_cord_paper/figures/Fig_5_line_plot_MN_count.pdf",
width = 4, height = 6,
plot = main
)
## Warning: Removed 53 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(
filename = "~/spinal_cord_paper/figures/Fig_5_ctrl_vs_poly_barplot_MN_count.pdf",
width = 3, height = 10,
plot = bar_cond
)
ggsave(
filename = "~/spinal_cord_paper/figures/Supp_Fig_5_ctrl_vs_poly_boxplot_MN_count.pdf",
width = 6, height = 5,
plot = box_day
)
pdf("~/spinal_cord_paper/figures/Fig_5_barplot_split_MN_count.pdf", width = 3, height = 10)
bar_axis
boxpl_axis
boxpl_0_axis
dev.off()
## svg
## 2
ggsave(
filename = "~/spinal_cord_paper/figures/Fig_5_boxplot_split_MN_count.pdf",
width = 7, height = 10,
plot = box_axis
)
# Date and time of Rendering
Sys.time()
## [1] "2025-04-16 14:38:05 CEST"
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: FlexiBLAS OPENBLAS; LAPACK version 3.11.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Zurich
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.3.0 ggpubr_0.6.0 ggplot2_3.5.1 tidyr_1.3.1
## [5] magrittr_2.0.3 stringr_1.5.1 tibble_3.2.1 dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.9 utf8_1.2.4 generics_0.1.3 rstatix_0.7.2
## [5] stringi_1.8.4 lattice_0.22-6 digest_0.6.37 evaluate_1.0.1
## [9] grid_4.4.1 fastmap_1.2.0 Matrix_1.7-0 jsonlite_1.8.9
## [13] backports_1.5.0 mgcv_1.9-1 purrr_1.0.2 fansi_1.0.6
## [17] scales_1.3.0 textshaping_0.4.0 jquerylib_0.1.4 abind_1.4-8
## [21] cli_3.6.3 rlang_1.1.4 munsell_0.5.1 splines_4.4.1
## [25] withr_3.0.2 cachem_1.1.0 yaml_2.3.10 tools_4.4.1
## [29] ggsignif_0.6.4 colorspace_2.1-1 broom_1.0.6 vctrs_0.6.5
## [33] R6_2.5.1 lifecycle_1.0.4 car_3.1-2 ragg_1.3.2
## [37] pkgconfig_2.0.3 pillar_1.9.0 bslib_0.8.0 gtable_0.3.6
## [41] glue_1.8.0 systemfonts_1.1.0 xfun_0.49 tidyselect_1.2.1
## [45] rstudioapi_0.16.0 knitr_1.49 farver_2.1.2 htmltools_0.5.8.1
## [49] nlme_3.1-165 labeling_0.4.3 rmarkdown_2.29 carData_3.0-5
## [53] compiler_4.4.1